Time evolution of projected entangled pair states in the single layer picture 
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We propose an efficient algorithm for simulating quantum many-body systems in two spatial 
dimensions using projected entangled pair states. This is done by approximating the environment, 
arising in the context of updating tensors in the process of time evolution, using a single-layered 
tensor network structure. This significantly reduces the computational costs and allows simulations 
in a larger submanifold of the Hilbert space as bounded by the bond dimension of the tensor 
network. We present numerical evidence for stability of the method on an antiferromagnetic isotropic 
Heisenberg model where good agreement is found with the available accurate results. 



I. INTRODUCTION 



Tensor network formalisms have been very successful in 
describing quantum many-body systems, and their the- 
oretical study is expected to play a crucial role for the 
understanding of strongly correlated phenomena in con- 
densed matter physics. Despite the huge Hilbert space 
associated to the many-body system, scaling exponen- 
tially with the number of particles, the introduction of 
the density matrix renormalization group [3 12] and ma- 
trix product states algorithms to simulate ground states 
of one dimensional quantum many-body systems [3l |4] 
have given strong evidence for the fact that physically 
interesting states are confined to reside in a small sub- 
manifold of the full Hilbert space. This observation was 
later proven in the context of quantum information the- 
ory in terms of entanglement properties of ground states. 
More specifically, it was shown that ground states of one 
dimensional systems whose Haniiltonian is gapped are 
only weakly entangled (they obey the area law) and can 
as such be faithfully and efficiently simulated in terms 
of matrix product states [S]. The main point of the 
tensor network formalism is that many-body quantum 
states are described in terms of local tensors (or matri- 
ces) where the number of parameters scales polynomially 
with the system size, and this in turn makes the classical 
simulations tractable. The matrix product state formal- 
ism was later generalized to two spatial dimensions by 
the introduction of 2 different methods, the multiscale 
entanglement renormalization ansatz (MERA) (BHEj and 
the projected entangled pair states (PEPS) O [10] both 
of which have later also been extended to fermionic sys- 
tems pTHTT] . Tensor network methods have a great po- 
tential in describing two dimensional quantum systems 
as they do not suffer from the notorious "sign problem" 
which makes frustrated spin systems and fermionic spin 
systems essentially untractable by quantum monte carlo 
methods. However, both MERA and PEPS suffer from 
relatively high computational complexity. The perfor- 
mance of PEPS algorithms is furthermore hindered by 
the instabilities that arise due to the lack of a normal 
form of the PEPS structure; this has to be opposed to 
the optimization using matrix product states in one di- 
mension where such a normal form makes the environ- 



ment essentially disappear from the local optimization 
procedure. 

The environment plays an important role in the con- 
text of both density matrix renormalization group and 
tensor network methods as it is responsible for a faithful 
projection of quantum states to a bounded submanifold 
of the Hilbert space, technically known as truncation. 
Still, the complexity of the PEPS algorithm scales as 
0[D^^) for open boundary conditions which essentially 
restricts the computations to small bond dimensions of 
D < 5 for finite size PEPS algorithms (but D < 8 for 
infinite PEPS algorithms) which is much smaller than in 
the one dimensional systems where one can reach bond 
dimensions of a few thousands. Note however that that 
much smaller bond dimensions should give already rea- 
sonable results; this follows from the monogamy proper- 
ties of entanglement. Another crucial difficulty in sim- 
ulating two dimensional many-body systems in terms of 
PEPS is however not only the computational complex- 
ity itself but also instabilities which occur due to the 
lack of the normal PEPS structure as is the case in one- 
dimensional systems. It turns out that the objects in the 
PEPS algorithm become less and less conditioned with 
increasing bond dimension and cutting ill-defined com- 
ponents immediately induces effective reduction of the 
bond dimension. This requires drastic changes to the 
PEPS simulation techniques, not so much in terms of 
complexity but rather in terms of stability, calling for 
the elimination of the double layer structure which is the 
root of both stability and complexity issues. The first 
step in this direction was made by not calculating the 
environment at all but letting it evolve in the process of 
imaginary time evolution |18) . Such an approach allows 
to achieve very large bond dimensions jl9| but neverthe- 
less seems to require a good initial approximation and 
works best for translation invariant PEPS states. 

In this paper we provide answers to both questions 
concerning the complexity and stability of the PEPS al- 
gorithms. We propose a method to simulate the time evo- 
lution of PEPS state using an approximate effective en- 
vironment, avoiding manipulations with the double layer 
tensor network. The approach is justified by the fact 
that the effective environment indeed exists and repro- 
duces the effect of the full environment on the system 
exactly. Similarly as in the density matrix renormaliza- 
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tion group method, the environment is approximated by 
a system of particles of an effective physical dimension 
which faithfully describe the effect of the environment to 
the system. The approach allows for high bond dimen- 
sions Z? = 10 or more, although the calculation of energy 
and other observables still requires the full calculation 
with the double layer tensor network structure. Alter- 
natively, quantum monte carlo sampling can be used to 
calculate the energy and expectation values of tensor net- 
works with a large bond dimension [l9Vl21j . 

II. METHOD 

A projected entangled pair state (PEPS) on a rectan- 
gular m X n lattice of qubits is parametrized in terms 
of local tensors ^I^^^I^^j for sites € {!,..., m} x 

{1, . . . , n} with local physical configurations Sij as 

E Tr(A[i'il^i'^...A["'"l^-")|si,i...s™,„) 

a) 

where tensors are contracted along the corresponding 
horizontal and vertical bonds between neighboring sites. 
The bonds connecting tensors across horizontal and ver- 
tical bonds between the sites are of dimension D such 
that aI^'^I"-'^ e C^xDxDxD^ rpj^g symbol Tr denotes 
the "tensorial trace" and implies contraction along the 
boundaries of the square. We shall assume open bound- 
ary conditions such that Tr will only refer to a map from 
a high-rank tensorial object (of unit size) resulting from 
the contraction of the square, to a scalar. 

Let us consider a bipartite splitting of the m x n sys- 
tem of qubits to a part consisting of a contiguous block of 
M qubits and a part containing the remaining [mn — M) 
qubits. The first part we shall call the subsystem S and 
the second part the environment E. If the subsystem S 
is subject to a local transformation resulting from e.g. a 
Suzuki- Trotter decomposition of the evolution operator, 
the internal bond dimension between the qubits in the 
subsystem S will increase. In order to keep the tensor 
network description |l]) manageable, the tensors A'*'-*' 
for (i, j) £ S must be truncated such that no bond di- 
mension exceeds the chosen bond dimension D. This is 
where the environment E comes into play due to the en- 
tanglement with the subsystem S. 

In this section we shall first show how the environ- 
ment can be efficiently approximated by an effective en- 
vironment which has approximately the same effect to 
the subsystem S. In the following we shall use the effec- 
tive environment to manipulate the subsystem S to lower 
the internal bond dimension. This is the core element in 
the temporal evolution of PEPS states where the evo- 
lution operator is decomposed using the Suzuki- Trotter 
decomposition into a product of local gates. 

For concreteness (and without lose of generality) let us 
consider a bipartition of the system where the subsystem 
S only contains two neighboring qubits in the horizontal 



direction, let us call them fj, = (J, J) and /i' = (/, J -f 1), 
which are for simplicity assumed not to be lying on the 
system boundary. If the subsystem S had been subject to 
a nontrivial local transformation, then the bond dimen- 
sion between sites /i and /i' has increased which calls for 
a truncation with the help of the environment consisting 
of all other sites. Let us contract the tensor network of 
the environment tensors {A^'''''" ,1/ £ E} which results in a 

joint environment tensor e}'^''^ '""^ where se = (si/, v £ E) 
denotes the physical (many-body) configuration of the 
environment sites. The PEPS state ([T]) is now rewritten 
to a compact form 

1^) = ^ Tr(A[^l''A[^'l'^'^['''^'l''''=)|s,s')|sE) (2) 

Despite exponentially large physical dimension of the en- 
vironment, Se G {1,2,..., 2™"~^}, it is only connected 
to the system S through six virtual bonds of dimension 
-D, adding up to a polynomially scaling bond dimension 
. Let us now decompose the state of the whole system 
to two parts by introducing an over complete set of states 
spanning the subsystem 5*, 

s,s' ,c 

whereas the environment is written as a superposition of 
configuration states in the environment as \'4'fiuu' dd' r)) ^ 

^tuJZM^)- The PEPS state ^ now takes a sim- 
ple form 

I*) =El^/)l^f) ^itl^ j = (WrfdV). (4) 

Due to the entanglement between the subsystem 5" and 
the environment, the former is in a mixed state given by 
the reduced density matrix ps = tr£;|4') (^Pj which reads 
(up to a normalization factor) 

pscx^(^f|^f)|^f)(Vf| (5) 

If one is to apply a local transformation to the subsystem 
S, such as a Trotter gate, the only relevant quantity to 
consider is the reduced density operator ps which in turn 
only depends on the environment through inner products 
{'fpj^\tpf) and not the state of the environment itself. This 
leads to the conclusion that, for a fixed set of basis states 
for the system S, there exists an effective environment 
of physical dimension which exactly reproduces these 
inner products and is given by 

\^f)- E^^"''"I^~e) (6) 

se = 1 

where E^^'^ ' is obtained by an orthogonal factorization 
e.g. E = QR where [E](,^),(i,,„„,dd') = E^^ruJZ' and 
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FIG. 1. Single layer contraction of PEPS structure (a): con- 
traction of two rows over vertical bonds (b) is followed by the 
truncation of the physical bond (c) resulting in (d). 



^\irul'dd') = [R-](5E),(ir-««'dd'); the unitary matrix Q is 
thus irrelevant. Such an effective environment suggests 
that the approximate contraction of the tensor network 
should be done already on the level of quantum states 
i.e. in the single layer picture, by truncating not only 
the virtual but also the physical degrees of freedom of 
lesser importance. 



A. Single layer contraction 

Let us show how the effective environment E^^'^ ' can 
be determined efficiently as depicted on Fig. [l] The con- 
traction takes place on the single layer with a bond di- 
mension D (as compared to for a double layer struc- 
ture) and results in an effective environment consisting 
of particles with a chosen effective physical dimension d 
connected by virtual bonds of dimension D. The tensor 
network is contracted row by row (or column by column) 
starting from above and from below such that the I-th 
row is surrounded by a single row of effective particles on 
both sides. 

The first step in the contraction scheme is well known 
|22j in the framework of tensor networks and involves 
contraction of two rows into a single row with an en- 
larged horizontal dimension and essentially squared phys- 
ical bond dimension. In order to make the method effi- 
cient, the resulting row must be truncated to a row with 
a bounded bond dimension D and a bounded physical 
dimension d. Let us consider a single layer contraction 
where rows are merged row-by-row starting from both 
upper and lower-most row. The result of contracting 
over vertical virtual bonds connecting two rows results 
in a matrix product state 

m = Y. tr(R™-^"^ • • ■ R["i^'"'''")isi,^'i, ■.■,s^,v^) 

(7) 

with a three- fold external bond dimension {sj,Vj) = 
{sj,s'j,Vj) where Sj and s'j are physical bonds at sites 
and (2,j), respectively, and Vj is the vertical bond 
connecting site (2,j) to site (3,j). The matrix product 
state I?! with large matrices R^'^-i''^ g qDDxdd 



easily be approximated by a matrix product state |4>) 
1$) - ^ tr(RW^i''^ • • • R["l^™'''")|si,z;i, . . . , 

with smaller matrices Rl-'l-i^^ e qDxd g^^j^ ^^ja^ the 
Euclid distance |||$) — |$)||2 is minimal. However, the 
physical bond dimension remains dd instead of the ini- 
tial d and continuing the procedure would result in an 
exponentially growing physical bond dimension. There- 
fore, the physical bond dimension must be truncated 
as well as depicted on Fig. If the matrix product 
state (|8| is written in an equilibrated form f3 , i.e. such 
that any given site in a row is connected to unitary en- 
vironments on both sides with the Schmidt coefhcients 
explicitly given on the corresponding bonds, a very good 
approximation to the optimal splitting (which is a quar- 
tic problem) can be found by truncating the physical di- 
mension at each site j independently^ by finding matrices 
Bbl*" for which ^2^^^'^'" ® BI-^I™ * best approximates 
RI-'I-^ (8) R[^1-" * according to the Frobenius norm. 
Such matrices are easily found from the singular value 
decomposition R, 

I^'" = T.UI^^'^^VI§, (9) 

S 

as B[J>~^ = If all singular vectors of U were 

retained where [U](/rt)),s — uj^^^^, such transformation 
would be exact while a good approximation to the envi- 
ronment is obtained taking d leading singular vectors of 
U in the singular value decomposition In the end, 
the two rows are described as a single matrix product 
state with a physical dimension d and a vertical external 
bond dimension D which guarantees bounded matrices 
and thus makes the algorithm efficient. While such ap- 
proximation is not strictly optimal, it nevertheless pro- 
vides a very good approximation in practice. The idea 
of truncating the physical degrees of freedom is not new 
but it is intrinsic to e.g. MERA [6 and also appears in 
other renormalization algorithms |23) . 

The procedure from the previous paragraph is repeated 
for all rows i < I and i > I starting from the upper 
{i — 1) and the lower [i — m) boundary row, respec- 
tively, such that the structure depicted on Fig. [T]i is ob- 
tained. Eventually, the same procedure is applied in the 
horizontal direction such that the two sites of interest are 
surrounded by ten (or less for boundary sites) effective 
environmental sites as shown on Fig. [2^. An effective 
PEPS structure obtained by this procedure (Fig. [2^) can 
be further simplified by absorbing corner sites to their 
neighbors (Fig. [2|d) which is done in a trivial way fol- 
lowed by an exact reduction of the effective physical bond 
(Fig#). 

This way an effective environment E^^'^ ' ^ is obtained 
which is in a straight-forward way related to the effective 




FIG. 2. Effective PEPS after the single layer contraction (a). 
Effective environmental sites on corner are absorbed in neigh- 
boring sites (b). The environment is approximated by two 
disconnected parts (c). 



double-layer norm operator A/'''*''' ' as 



_yy[A',A''l _ ^ £'[m,A'']se* 



>E 



[m,m']se 



(10) 



where its matrix elements represent the inner products 
Aflj''^ ] ^ i^^E^^Ej fjgfine;;! previously. In the context 
of the usual time evolution, the effective operator A/" is 
used to determine a new set of tensor {Al'^l,^!^ 1} after 
a Trotter step has been applied on sites (/i,/i'). 

Similarly to the double layer contraction, the sin- 
gle layer contraction becomes exact for sufficiently large 
truncation parameters D and d whereas good approxima- 
tions can be obtained already with d — D = D. On the 
other hand, this way of calculating the environment is in 
many ways advantageous to the conventional contraction 
of the double layer tensor network. The first advantage is 
the computational cost. While the complexity of the con- 
ventional double layer contraction scales as 0{D^^), the 
costs of the single layer truncation only scale as 0{D''). 
The second advantage is that by doing the single layer 
truncation, good estimates can be obtained of how to 
perform gauge transformations on the original environ- 
ment sites connected to the subsystem S, that the ef- 
fective environment norm M becomes better conditioned 
[25j . which is of crucial importance for the stability of 
the algorithm. 

Let us now consider a Trotter gate T^^'^ 1 acting on 
two neighboring sites (/j,, and find a matrix product 
state 5* defined as 



I*) 



E 



Tr M^"' A 



(11) 



which best approximates T^^^ [vj/) such that the bond 
dimension between the tensors on sites (/i, /i') is upper- 
bounded by D. In the conventional update scheme, 
the tensors {A^^l,^!^ 1} are obtained by solving a multi- 
quadratic optimization problem in an iterative way which 
involves solving a linear system of equations N^A^^^^ = 
where the iv]^^' is an effective norm operator for the site 
/i, obtained by contracting M^^'^ 1 and tensors at the site 
n' (and similarly for site ^') . While the gauge trans- 
formations allow us to make iV^jg equal to the identity in 
the case of matrix product states, this is not possible in 
the case of PEPS states. Furthermore, the linear system 
of equations in consideration is typically ill-conditioned 




FIG. 3. A PEPS state with an effective environment (a) (iden- 
tical to Fig. |2|d) can be understood as a matrix product state 
with periodic boundaries (b). 



due to emergence of very small eigenvalues in the effec- 
tive norm operator N^s leading to instabilities, especially 
for large bond dimensions I? = 3,4, which become more 
and more pronounced in the process of simulation as we 
shall show later. In the usual PEPS time evolution, this 
problem is evaded by the projection to the "well" de- 
fined subspace spanned by the singular vectors of N^g 
with respect to a cut-off parameter e determining the 
ratio between the smallest kept and the maximal singu- 
lar value. If the parameter is set too high, the effective 
bond dimension is largely reduced whereas in the case of 
too low setting some ill-conditioned components are also 
retained resulting in instabilities. While one could use 
^[m.m ] obtained by the single layer contraction to calcu- 
late the effective norm operator (10 1 and then use the 



conventional update scheme, that would still not solve 
the stability issues. 



B. Environment splitting 

In order to eliminate the stability problems, we will 
show how to avoid calculations with the double layer 
tensors such as N defined (10 1. The core of the prob- 



lems is that the environment appears as a cyclic ma- 
trix product operator (see Fig. ^p) and as such does not 
permit the standard way of finding the optimal tensors 

using the singular value decomposition. 
We however know empirically, that the environment of 
two sites in sufficiently large lattices can be fairly well 
approximated by a product of two separate environments 
(Fig. [2]:) which is the idea we will pursue in the following. 

Let us rewrite the PEPS state ([2]) as a matrix product 
state with periodic boundaries (see Fig. |3]) as 

l^-) =^tr(L^-A['^l^A['^'l^'R^-)|sL,s,s',Si^) (12) 



where [A[''1^]( 



\(lud)r - ^Irud^ [J^^'^ ' Uiudr) - Arud 

whereas L" and correspond to the contraction of ten- 
sors belonging to the left and right three sites on Fig.[2]D. 

There is no way known to us how to found the optimal 
matrices aI'^J" and A^'^ 1" exactly without employing the 
sweeping optimization mechanism described in the previ- 
ous section. However, there are several ways to split the 
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environment approximately, assuming that the internal 
correlations between the two parts of the environment 
decay sufficiently fast. The first possibility which gives 
remarkably good results is to do the singular value de- 
composition of both parts of the environment, 



[L] y[i]T/[i]^ 



(13) 



and then taking the left approximate environment as 



if 



s ^ Ls 



(14) 



The right part of the environment is transformed in a 
similar way. In practice, the singular value is done sep- 
arately for values of the internal environment bond 7, 
followed by the singular value decomposition of the con- 
catenated and weighted right singular vectors which is 
numerically favorable. 

The approach in the previous paragraph is rather ex- 
pensive in our case due to the large physical dimension 
of the environment sites [d?). For that reason, we will 
pursue in an even simpler way by simply self-contracting 
the internal environment bond for each part of the envi- 
ronment 



(15) 



The result is a single matrix product state with open 
boundary conditions for which the optimal matrices 
and At'^ I'* are determined exactly by the singular value 
decomposition. From numerical tests we observe that 
this approach is only slightly less accurate than the one 
presented earlier. The total cost of this step scales as 
O(-D^), all coming from the SVD, but it is in practice neg- 
ligible for relatively small dimensions D ^ 10 when com- 
pared to the single layer contraction part of the method 
involving fairly many steps scaling as 0(1?^). 



III. RESULTS 

To illustrate the validity of the method we consider an 
antiferromagnetic Heisenberg model on a square lattice 



(16) 



for which the ground state properties are accurately de- 
scribed by the stochastic series expansion (SSE) quan- 
tum nionte carlo method |24| . We perform the imaginary 
time evolution \^{(3)) = e~^^\^Q) using the second or- 
der Suzuki- Trotter expansion where two-site local gates 
are applied to the PEPS state, followed by the truncation 
of the corresponding tensors. 

First we compare the proposed method to the usual 
imaginary time evolution of PEPS [10 with a bond di- 
mension _D = 2 on a system of 6 x 6 qubits. In the 
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FIG. 4. Imaginary time evolution of a random initial state 
with bond dimension D = 2 under Hamiltonian ( 16 1 on a 



6x6 system using the usual PEPS time evolution algorithm 
(ITE) and the single layer time PEPS algorithm (SLTE). The 
Trotter time step was set to 10""^. The cutoff in the linear 
problem e appearing in ITE is designated in the legends. The 
D = d for the SLTE is given in the legends whereas D = 64 for 
the ITE. The exact ground state energy equals Eq = —86.9. 



single layer method we choose D = d ~ 4,6 whereas 
the cut-off parameter in the usual PEPS time evolution 
is set to Dcut = 64. In the latter, ill-conditioned lin- 
ear systems require another cut-off parameter which is 
chosen as e = 10~^°,10~^. The results in Fig. [4] con- 
firm that the usual time evolution results in instabilities 
due to the ill-conditioned linear problems which are being 
solved in the simulation. Increasing the cut-off parameter 
£ from 10"^'' to 10~^ suppresses the non-physical solu- 
tions and pushes the simulation forward, although the 
accuracy of simulation steps is reduced. Since the deter- 
mination of the cut-off parameter is heuristic procedure 
and inevitably results in either cutting relevant degrees 
of freedom or keeping non-physical ones, the simulation 
eventually becomes unstable. This phenomenon is even 
more pronounced with larger bond dimensions where the 
linear problems are of larger dimension and it is even 
more difficult to make a sensible compromise for s. From 
the technical point of view, the PEPS tensors are always 
rescaled such that their 2-norm is the same for all ten- 
sors. Increasing the cutoff parameter D to 128 did not 
help significantly for e = 10"^*^. 

The single layer time evolution, on the other hand, pro- 
duces no instabilities for arbitrarily long times, although 
the results do oscillate slightly (not noticeable on the fig- 
ure). Choosing a larger effective bond dimensions D and 
d makes the convergence faster but nevertheless results in 
an comparably good final state. We note, however, that 
the single layer time evolution involves approximations 
and slightly more accurate results can be achieved (note 
a few points for ITE, e = 10~^) using the usual time evo- 
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FIG. 5. Absolute error of the energy per site for the antifer- 
romagnetic Heisenberg model ( 16 1 as a function of iteration 



step and lattice sizes L x L as given in the legend. The results 
for D = 6 were obtained by sampling. 



lution. The time needed to obtain the results by the usual 
time evolution is by a factor of hundred larger than the 
single layer time evolution. Both approaches were done 
using the same equally (non-)optimized code and start- 
ing from the same random initial state. For smaller bond 
dimensions, the single layer approach can thus be used to 
quickly obtain a very good initial approximation which 
could be further refined by some double layer technique 
such as minimizing the energy by sweeping [9]. 

Secondly, we test the method on larger systems and 
larger bond dimensions. For a given bond dimension, 
the simulation was running as long as the energy decay 
rate was sufficiently high and the results were used as the 
initial state for simulations with a larger bond dimension 
D + 1. The growth from D to D + 1 is intrinsic to the 
time evolution and does not require any zero-padding. 

In Fig. [5] we present the results for the absolute error 
of the energy per site for the antiferromagnetic Heisen- 
berg model ([16]) compared to the results obtained by the 
stochastic series expansion which we take as exact. We 
consider three system sizes and bond dimensions up to 
D — 6. The single layer truncation parameters were in 
all cases chosen as D ^ d = D. Note that this model 
is critical and is among fairly difficult models to sim- 
ulate using tensor network methods. From the results 
for D < 5 where the energy is calculated by the approxi- 
mate contraction of the PEPS tensor network, it is clearly 
visible that the energy decreases monotonically until the 
plateau is reached which justifies the single layer contrac- 
tion scheme and the approximate environment splitting. 

We can easily simulate PEPS systems of 10 x 10 sites 
with bond dimensions _D = 10 or even more, however 
the extraction of expectation values such as the energy is 
nontrivial unless we also transform the hamiltonian itself 
by the isometrics generated in the single layer contrac- 



tion procedure as known in the context of the DMRG. For 
purposes of this manuscript we rely on the double layer 
contraction scheme to calculate the energy which, much 
more than in the translation invariant infinite PEPS algo- 
rithm, is computationally very costly and at present only 
sensible for bond dimensions D ^ b. For that reason, 
the energies of the PEPS with a large bond dimension 
D = Q are calculated by sampling using the Metropolis 
algorithm following the fact PHI HI] that 



E 



where = (17) 



Metropolis algorithm also allows us to sample not j^*) 
but rather P\s^=o\'^) where P\s^=() is a projection op- 
erator to the zero total spin subsector. The benefit of 
the Metropolis sampling is again the single layer picture 
of the problem which avoids doubling the virtual bond 
dimension and thus reduces the contraction cost. The 
downside however is that one has to carefully choose the 
Metropolis updates otherwise the variance decays rela- 
tively slowly as visible from Fig. [5] where the variance is 
still large. However, the results obtained by sampling for 
D = 6 are consistent with the results for Z? = 5 obtained 
by the approximate contraction of the tensor network. 



IV. CONCLUSION 

We have presented an approximate method to simulate 
quantum many-body systems on finite two-dimensional 
lattices using the projected entangled pair states (PEPS) 
where the contractions are only done on the level of quan- 
tum states (single layer tensor network). This contrasts 
with the conventional PEPS simulation scheme where the 
simulations involve contractions on the level of expecta- 
tion values given by a double layer tensor network with a 
squared bond dimension, leading to a high computational 
cost. The single layer approach eliminates the stability is- 
sues present in both the time evolution of PEPS and the 
variational PEPS algorithm to simulate ground states, 
which in both cases originate from the ill-conditioned 
double layer effective norm operator. Unlike the usual 
time evolution which is exact for sufficiently large cut- 
off parameters and arbitrary precision arithmetics, the 
single-layer time evolution is approximate due to its sim- 
plified treatment of the environment. We compared the 
single layer time evolution to the usual time evolution 
for PEPS and observed comparable results in accuracy 
whereas the stability of the simulation is better in the sin- 
gle layer approach. We tested the method for larger bond 
dimensions on an isotropic antiferromagnetic Heisenberg 
model on a square lattice where we again observed mono- 
tonic decrease of the energy for bond dimensions D < 5 
where it can be calculated by contracting the PEPS ap- 
proximately. We presented sampled results for D — 6 
which are consistent with the results for lower bond di- 
mensions. While we only presented the results for a con- 
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ceptually simple model, the method can easily be applied 
to any spin or fermionic system in two spatial dimensions. 
The extension of the single layer contraction technique to 
the latter case is straight-forward where all arising sign 
factors can be absorbed locally due to the well defined 
parity of tensors of the fermionic PEPS |TT| . 



tational results presented have been achieved using the 
Vienna Scientific Cluster (VSC). 
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